function F = cal_chg_bartik(param, const, eqm)
%% Choose optimal quality after Bartik shock
        
    % D^B(q,E=1)
    rv_B = (eqm.D_H/param.f_vH).^(1/(param.beta_v-1))./( (eqm.D_H/param.f_vH).^(1/(param.beta_v-1))+(param.e^param.sigma*1.05*const.D_F/param.f_vF).^(1/(param.beta_v-1)) );
    D1_B = rv_B.*eqm.D_H + (1-rv_B)*param.e^param.sigma*1.05.*const.D_F;        

    % f_v^B(q,E=1)
    f_v1_B = rv_B.^param.beta_v*param.f_vH + (1-rv_B).^(param.beta_v)*param.f_vF;

    % C_v^B(q,E=1)
    C_v1_B = f_v1_B.^(-1/param.beta_v).*(1./(param.sigma*param.w)).^(1/param.beta_v);

    % C_x^B(q,E=1,I)
    C_x11_B = ((param.sigma/(param.sigma-1)*param.w.^(1-param.alpha_m-param.alpha_s)).^(1-param.sigma) ...
               .*eqm.C_m1.^param.alpha_m.*C_v1_B).^const.gamma;
    C_x10_B = ((param.sigma/(param.sigma-1)*param.w.^(1-param.alpha_m-param.alpha_s)).^(1-param.sigma) ...
               .*const.C_m0.^param.alpha_m.*C_v1_B).^const.gamma;
      
    % Pi^B(q,E=1,I)
    Pi11_B = D1_B.^const.gamma.*(eqm.c1.^param.alpha_m*param.P_s^param.alpha_s).^((1-param.sigma)*const.gamma).*C_x11_B;
    Pi10_B = D1_B.^const.gamma.*(eqm.c0.^param.alpha_m*param.P_s^param.alpha_s).^((1-param.sigma)*const.gamma).*C_x10_B;

    % Import Cutoff and Probability
    lnZQ = const.k0*const.lnzQ-log(param.sigma*const.gamma);
    ImportCutoffQ1_B = (lnZQ+log(kron(ones(param.Nf,1),Pi11_B)-kron(ones(param.Nf,1),Pi10_B))-param.mu_I)/param.sigma_I;
    ImportProbQ1_B = normcdf(ImportCutoffQ1_B);  

    % E^B[pi(q,omega,E=1)]
    EPiQ1_B = exp(lnZQ).*(kron(ones(param.Nf,1),Pi10_B) + kron(ones(param.Nf,1),Pi11_B-Pi10_B).*ImportProbQ1_B) ...
        - exp(param.mu_I+1/2*param.sigma_I^2).*normcdf(ImportCutoffQ1_B-param.sigma_I);

    % NetEPiQ
    NetEPiQ_B = EPiQ1_B-eqm.EPiQ0;
    if min(min(NetEPiQ_B))<0
        disp('Adjustment: EPiQ1_B-EPiQ0 < 0')
        NetEPiQ_B(NetEPiQ_B<0) = 0;
    end 

    % Export Cutoff and Probability
    ExportCutoffQ_B = (log(NetEPiQ_B)-param.mu_E)./param.sigma_E;
    ExportProbQ_B = normcdf(ExportCutoffQ_B);     

   % Choose new quality (grid search)
    [qB_index, ~] = grid_search(param, eqm.EPiQ0, EPiQ1_B, ExportProbQ_B, ExportCutoffQ_B);        


      
%% Calculate average wage response

    % Quality level
    q_level = const.Qgrid(eqm.q_index)'; 
    qB_level = const.Qgrid(qB_index)';
    
    % Firm wage
    wageO = const.wage_grid(eqm.q_index)';
    wageB = const.wage_grid(qB_index)';   
    
    % Average quality and wage response  
    weight = eqm.density_exp;
    avg_quality_response = sum(qB_level./q_level.*weight)/sum(weight);    
    avg_wage_response = sum(wageB./wageO.*weight)/sum(weight);
 
 
%% Changes in average wage of suppliers   
  
    avg_wgt_lnwageS_q = eqm.theta_m.*eqm.c_H.^(param.sigma-1)./eqm.V.*sum(const.phi_y.*const.phi_v.*kron(ones(param.Q,1), (log(param.w)+const.lnwage_grid).*eqm.Psig),2)'; 
    avg_wgt_lnwageS = avg_wgt_lnwageS_q(eqm.q_index)';
    avg_wgt_lnwageSB = avg_wgt_lnwageS_q(qB_index)';
    avg_wageS_response = sum((avg_wgt_lnwageSB - avg_wgt_lnwageS).*weight)./sum(weight);
   
    
%% Export intensity

    %rB = (eqm.D_H/param.f_vH).^(1/(param.beta_v-1))./((eqm.D_H/param.f_vH).^(1/(param.beta_v-1))+(param.e^param.sigma.*eqm.D_F.*param.shock/param.f_vF).^(1/(param.beta_v-1)));
    %DB = rB.*eqm.D_H + (1-rB).*param.e^param.sigma.*eqm.D_F.*param.shock;     
    %export_intensity = 1-rB.*eqm.D_H./DB;        


%% Return

    %bartik.PiB                  = PiB;
    bartik.qB_index             = qB_index;
    bartik.avg_quality_response = avg_quality_response;
    bartik.avg_wage_response    = avg_wage_response;
    bartik.avg_wageS_response   = avg_wageS_response;
    %bartik.export_intensity     = export_intensity;
        
    F = bartik;

    
end